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ABSTRACT 

This study investigates the evolution of Rayleigh- Taylor (R-T) instabilities in 
Type la supernova remnants that are associated with a low adiabatic index 7, where 
7 < 5/3, which reflects the expected change in the supernova shock structure as a result 
of cosmic-ray particle acceleration. Extreme cases, such as the case with the maximum 
compression ratio that corresponds to 7 = 1.1, are examined. As 7 decreases, the shock 
compression ratio rises, and an increasingly narrow intershock region with a more pro- 
nounced initial mixture of R-T unstable gas is produced. Consequently, the remnant 
outline may be perturbed by small-amplitude, small-wavelength bumps. However, as 
the instability decays over time, the extent of convective mixing in terms of the ratio 
of the radius of the R-T fingers to the blast wave does not strongly depend on the 
value of 7 for 7 ^ 1.2. As a result of the age of the remnant, the unstable gas cannot 
extend sufficiently far to form metal-enriched filaments of ejecta material close to the 
periphery of Tycho's supernova remnant. The consistency of the dynamic properties of 
Tycho's remnant with the adiabatic model 7 = 5/3 reveals that the injection of cosmic 
rays is too weak to alter the shock structure. Even with very efficient acceleration of 
cosmic rays at the shock, significantly enhanced mixing is not expected in Type la 
supernova remnants. 
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1 INTRODUCTION 

Convective or Rayleigh- Taylor (R-T) instabilities occur in 
the intershock region of young supernova remnants (SNRs), 
whose evolution is dominated by the structure of the ejecta 
of the supernova SN. 

After the SN explodes, its fast ejecta are continuously 
decelerated by ambient material (interstellar or circumstel- 
lar medium). Because the flow is supersonic, a blast wave 
propagates ahead of the spherical front of the expanding 
gas, or contact discontinuity (CD.), whose pattern speed 
is equal to the flow speed of the ejecta. As the ejecta gas 
gradually piles up at the CD., a high pressure is established 
and an internal shock arises, which moves inward relative to 
the expanding ejecta in a Lagrangian sense. An SNR thus 
comprises two shells of heated, shocked gas, bounded by a 
forward and reverse shock, separated by a CD. The reverse 
shock eventually turns back into the explosion center and 
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terminates the free expansion of the supernova ejecta. In 
the intershock region, since the gas experiences outwardly- 
directed effective gravity, owing to deceleration, the flow be- 
comes convectively unstable, such that the gas in the inner 
shocked shell protrudes into, and mixes with, the gas in the 
outer shell. Such an instability is only expected during an 
early SNR phase of expansion, when the SN ejecta bear a 
sufficiently sharp density distribution relative to the ambi- 
ent medium, enabling the convection to persist (Wang & 
Chevalier 2001 & 2002; hereafter WC01 & WC02). 

Rayleigh- Taylor instability has been considered to be 
effective for modifying the spherical morphology of an SNR. 
A commonly observed SNR with irregular remnant outlines 
is Cassiopeia A SNR (Cas A), a prototype of Type II SNRs, 
of which both optical (Fesen & Gunderson 1996) and X-ray 
(Hughes et al. 2000) observations reveal numerous bumps 
and protrusions on its roughly-spherical outline. SN 1986 J 
and SN 1993 J are other, extreme examples of core-collapse 
SNe (Bartel et al. 1991; Bietenholz, Bartel, & Rupen 2001). 
Such an irregular morphology has been explained by the 
strong intershock convection that spreads the ejecta mate- 
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rial out and beyond the nominal radius of the forward shock, 
but the extensively mixed structure, including the dispersal 
of 56 'Ni to high velocities, is now believed to be imprinted 
by the R-T instability that occurs during an asymmetric 
explosion (Kifonidis et al. 2000, 2003; Hungerford, Fryer, 
& Warren 2003), although processes which occur after the 
explosion such as the nickel bubble effect, in which the ra- 
dioactive progenitors of the Fe expand relative to their sur- 
roundings, explaining the large filling factor of the observed 
Fe (Woosley 1988; Li et al. 1993), contribute to the mixing. 
With respect to young Type la SNRs, X-ray observations of 
Tycho's remnant (SN1572) reveal heavy ejecta material and 
dense knots, which, in non-decelerating motion, reach close 
to, and protrude from, the eastern periphery of the remnant 
(Hughes, 1997; Hwang, Hughes, & Petre 1998). Radio obser- 
vations of Tycho's remnant (Reynoso et al. 1997) similarly 
revealed that the northeastern and the adjacent southeast- 
ern parts protrude from the circular outline; the southeast- 
ern part corresponds to the X-ray knots, while a regularly 
spaced structure immediately beneath the sharply marked 
forward shock front is present at the northeaster edge (Ve- 
lazquez et al. 1998). The large extension and non-decelerated 
motion of these structures are attributed to dense clumps in 
the diffuse supernova ejecta that expand into the remnant 
shell, rather than the global R-T convection that is caused 
by deceleration of the remnant (WC01). 

Although various hydrodynamic simulation models 
have been applied to study R-T instability in SNRs, most 
use a power-law density distribution of the SN ejecta, which 
yields a self-similar solution for the shock structure and 
applies only to the early expansion phase of core-collapse 
SNRs, when the inner flat component of the ejecta has 
not expanded so far as to change the intershock structure 
(Chevalier, Blondin, & Emmering et al. 1992; Jun & Nor- 
man 1996a; Blondin & Ellison 2001, hereafter BE01, WC02; 
Fraschetti, 2010) Moreover, these idealized models with an 
adiabatic index of 7 = 5/3, such as that of Chevalier et 
al. (1992), have revealed rather limited R-T mixing. Jun 
& Norman (1996a and 1996b) reached a similar conclusion 
after considering rapid cooling (Chevalier & Blondin 1995) 
in the shocked ejecta. Jun & Jones (1999) used accelerated 
cosmic-ray electrons as test particles to calculate the ra- 
dio synchrotron emission but emphasized the reverse shock. 
Since rapid cooling and efficient cosmic-ray particle acceler- 
ation (Ellison, Berezkho, & Baring 2000) both enhance the 
shock compression of the gas, BE01 adopted a low adiabatic 
index 7 < 5/3 to examine the instability, which exhibited 
the change expected because of particle acceleration. In this 
case, the mixing extends far enough to reach the forward 
shock front and considerably perturb the remnant outline. 
Recently, Fraschetti et al. (2010) employed 7 = 4/3 to rep- 
resent the relativistic cosmic-ray particles and Ferrand et 
al. (2010) applied a kinetic model of nonlinear diffusive par- 
ticle acceleration to the same hydrodynamic model as in 
Fraschetti et al.; they found that the development of insta- 
bility is not strongly enhanced, independently of whether 
the cosmic rays dominate the reverse shock, or the shocks 
receive a back-reaction of cosmic rays that evolves with time. 

Notably, the structure used in the above hydrodynamic 
models is not appropriate for Type la supernovae (SNe la), 
whose density falls exponentially with increasing velocity 
and is significantly different from that given by the power- 



law model (Dwarkadas & Chevalier 1998; Dwarkadas 2000; 
WC01). An exponential profile is formed probably because 
the subsonic explosion of SNe la steadily releases energy be- 
hind a burning front, unlike a pure shock acceleration in 
core-collapse SNe. The greatest dynamical contrast between 
the two models is that, because the exponential ejecta den- 
sity gradient decreases with expanding SNR, the convec- 
tion in Type la SNRs decays over time. As a result of the 
age of the remnant, the mixture of metal-enriched gas near 
the blast wave of Tycho's remnant cannot be attributable 
to an intershock convection. Based on the indication that 
ejecta clumps are broken-up components of the shell of an 
inflated Ni bubble in the center of supernova ejecta that is 
heated by radiation from the 56 Ni -> 56 Co -> 56 Fe decay 
sequence (Basko 1994, WC01), the author (Wang 2005 & 
2008) utilized radiative transport radiation-hydrodynamic 
simulations to examine the structure and expansion proper- 
ties of the inferred clumps in core-collapse and Type la SNe. 
For Type la SNe, when the density distributions produced 
by successful, W7-like SN explosion conditions are applied to 
the ejecta substrate, the properties of the inferred clumps is 
maximally compatible with those of the X-ray knots and fil- 
aments near the outline of Tycho's remnant (Wang 2008). In 
contrast, Warren et al. (2005) argued that, given the Chan- 
dra X-ray evidence of efficient cosmic-ray particle accelera- 
tion in Tycho's blast wave, the observed extended mixture 
of ejecta material could be explained by the R-T convection 
in the remnant shell, as was proposed by BE01. These in- 
vestigations cast doubt on whether the enhanced intershock 
instabilities, expected to result from cosmic-ray particle ac- 
celeration, are responsible for the mixing present at Tycho's 
periphery. 

This study examines the R-T instabilities in Type la 
SNRs using a low adiabatic index to approximate the ex- 
pected increase in shock compression ratio caused by cosmic- 
ray particle acceleration. This study differs from the previ- 
ous work of BE01 in that it adopts an exponential ejecta 
model to describe the general structure of Type la SNe, 
allowing the SNR evolutionary state to change with time. 
The rest of this paper is organized as follows. Section 2 in- 
troduces the effect of particle acceleration on SNRs. Section 

3 describes the computational setup and methods. Section 

4 then elucidates the evolutions of the R-T instabilities in 
one and two dimensions. Finally, section 5 and 6 presents 
discussions and draws conclusions. 



2 ACCELERATION OF COSMIC-RAY 

PARTICLES AND ASSOCIATED SHOCK 
COMPRESSION RATIO 

Cosmic rays (CR) are electrons and atomic nuclei, mostly 
protons, that move at nearly the speed of light. Their en- 
ergies can exceed 10 21 eV. Strong shock waves of supernova 
remnants are believed to be efficient accelerators and pri- 
mary sources of Galactic cosmic-ray particles, and to boost 
the energies of such particles to the knee of the CR spectrum, 
10 15 eV (Axford 1981). The signature of this process is the 
decay of pions, which are generated in collisions of acceler- 
ated protons with atom and molecules while the shocks pass 
through the surrounding interstellar medium (ISM). Elec- 
trons are also known to be accelerated to cosmic-ray ener- 
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gies in supernova remnants, via the scattering of low-energy 
photons, mostly the 2.7K cosmic background photons, which 
become TeV gamma rays through inverse Compton scatter- 
ing (Koyama et al. 1995; Mastichiadis & de Jager 1996, and 
Tanimori et al. 1998 for observations of SN 1006). To ex- 
plain the spectrum of SN 1006, Chevalier & Reynolds (1981) 
first suggested that synchrotron emission, rather than other 
thermal and non-thermal processes, was responsible for the 
high-energy photons in SNRs. The energies of CR that are 
observed on Earth and the strength of the interstellar mag- 
netic fields deduced from radio measurements support the 
prediction that the synchrotron radiation of CRs has X-ray 
wavelengths. Observations of X-ray synchrotron radiation 
from electrons ever since have shown that SNRs produce 
GeV emissions, although the synchrotron radiation in radio 
and optical wavelengths is not linked with the high-energy 
CRs (for a review see Reynolds 2008). 

Although growing evidence relates cosmic-ray ions to 
the thermal properties of the shock-heated X-ray-emitting 
gas (Decourchelle, Ellison, & Ballet 2000; Ellison, Slane, & 
Gaensler 2001; Ellison et al. 2010), evidence of ion/proton 
acceleration in SNRs remains circumstantial. Cascade show- 
ers of optical photons from pion decay in RX J 1713.7- 
3946, which appear to be similar to SN 1006, even reveal 
the presence of TeV 7-rays (Enmoto et al. 2002). The non- 
thermal X-ray emission from the shells of Galactic super- 
nova remnants, most notably SN 1006 and Cas A, and the 
dominant non-thermal component in a few newly-discovered 
SNRs, believed to be the synchrotron emission from elec- 
trons that is accelerated to ~ 100 TeV, represents evidence 
that cosmic rays are produced in SNR shocks (Petre, Hwang, 
& Allen 2001). However, most other remnants appear to 
be dominated by thermal emission from their shells, and 
disagreement exists about whether the 7-ray emissions in 
the direction of RX J 1713.7-3946 and others come from 
ions (Berezkho & Volk 2006) or electrons (Katz & Waxman 
2008). Reynolds & Keohane (1999) calculated the maximum 
energies of shock-accelerated electrons in 14 young shell- 
type SNRs and concluded that young remnants could not 
be responsible for the highest-energy Galactic cosmic rays. 
A TeV proton has a 10~ 10 times lower critical synchrotron 
emission frequency than a TeV electron; however, identify- 
ing synchrotron emission from cosmic-ray protons in the ra- 
dio wavelength is impossible, considering that a TeV proton 
produces a ~ 10 -3 times smaller flux than a TeV electron 
and many more GeV electrons than TeV protons contribute 
to the emission. Even if a spectral signature of proton/ion 
acceleration could be found, for which the best and only di- 
rect evidence relies on the gamma-ray emission that is not 
associated with inverse Compton emission from electrons 
while the indirect evidence comprises measurements of ab- 
normally low electron temperatures and curvature in the 
SNR spectrum (Glenn Allen 2011, private communication; 
Allen, Houck, & Sturner 2008), the evidence in support of 
the claim that isolated SNRs are the main accelerators of 
galactic CRs may be insufficient (Butt 2009). 

To generate high-energy cosmic rays, the acceleration 
must transform at least 10% of the total ejecta kinetic en- 
ergy into relativistic particles, and this process cannot oc- 
cur too soon (~ 100 yr) after the supernova explosion, as in 
that case, the cosmic rays would lose most of their energy 
(Blandford & Eichler 1987; Drury, Markiewicz, & Volk 1989; 



Dorfi 2000). When efficient particle acceleration occurs, rel- 
ativistic, superthermal particles can be created. The efficient 
production of superthermal particles is accompanied by de- 
creased postshock temperature. Consequently, the region be- 
tween the forward and reverse shocks of an SNR shrinks and 
becomes denser. This effect increases the compression ratio 
of the shock from a — 4, where 



for strong shocks, to a > 4; and so the effective adiabatic 
index is less than 5/3 for the SNR intershock region. The 
increased compression of the gas gives rise to a more com- 
pact intershock structure and is expected to cause larger 
amplitudes for the R-T instability. 

In forecasting the gamma-ray flux from pion decay, 
Chevalier (1983) considered two-fluid, self-similar solutions 
where the ejecta was a thermal gas with 7 = 5/3 and the 
ambient medium was a relativistic gas that represented CRs 
with 7 = 4/3. In the case in which relativistic particles at 
the shock front contribute 100% of the total pressure, the 
effective adiabatic index is 1.333, and the shock compres- 
sion ratio is similar to that in the late radiative phase of 
SNR expansion, a > 7. BE01 showed that the results of the 
power-law model with 7 = 1.1 (a = 21) are similar to the 
case of radiative cooling that was considered by Chevalier 

6 Blondin (1995) where the ejecta behaved as an ideal gas 
with gamma=l and the ambient medium was modeled with 

7 = 5/3. In kinetic simulations of diffusive shock acceler- 
ation (DSA), the cosmic-ray-modified shocks depend only 
weakly on the Mach number Mo of the shock, and the total 
compression ratio is less than 10 even for Mq ~ 100, be- 
cause the propagation and dissipation of Alfven waves up- 
stream reduces the CR acceleration and the precursor com- 
pression (Kang, Ryu, & Jones, 2010). Based on these cal- 
culations, this work considers a maximum compression that 
corresponds to 7 = 1.1 for a young SNR such as Tycho's 
SNR, although theoretically arbitrarily large compression 
ratios are allowed. 

Berezhko & Ellison (1999) demonstrated that the de- 
pendence of the compression ratio on the sonic and Alfven 
Mach number implies that the compression is most pro- 
nounced in SNRs with large shock speeds but relatively weak 
magnetic fields. However, for a shock with a high Mach num- 
ber, a ~ 4 can still hold if the injection of cosmic-ray parti- 
cles into the shock is weak. 



3 MODELS AND METHODS 

The two-dimensional hydrodynamic simulations performed 
in this investigation is based on the the two-dimensional 
finite-difference code ZEUS2D (Stone & Norman 1992), fol- 
lowing the methods adopted in the author's previous work 
WC01. An exponentially-declining ejecta and a constant am- 
bient interstellar medium is initialized on a spherical grid. 
The flow velocity of the ejecta is freely-expanding and the 
ambient ISM is at rest. The inner boundary is inflow and 
fixed in radius, while the grid is radially nonuniform and 
expanding, following the motion of the forward shock and 
allowing finer resolution in the intershock structure. A grid 
of 600 radial zones typically resolves the intershock region 
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into ~ 400 zones. Simulations usually begin as early as 10 s 
after the supernova explosion, such that the flow can become 
numerically saturated and fully developed into the nonlinear 
regime long before Tycho's present epoch. Simulations initi- 
ated earlier presented numerical difficulty because the higher 
density reduced the numerical time step determined by the 
Courant condition. The pressure distribution of the ejecta 
satisfies the adiabatic law p = up"' , where ft = 6.1 X 10 14 
(cgs units), if the thermal energy is provided mainly by the 
entropy change during the nuclear burning of C and/or O 
to 56 Ni and 7 = 4/3 (Wang 2005). In the ejecta-dominated 
phase, the Mach numbers of the forward and reverse shocks 
greatly exceed 1, and so the gas pressure used in the ejecta 
and ambient medium was trivial. 

The exponential density profile of the ejecta is charac- 
terized by an explosion mass M and energy E, described 

by 

Psn = Aexp(-v/v e ) t~ 3 , (2) 

where v — r/t is the flow velocity of the ejecta, while A is 
a constant, A = 6 3/2 M 5/2 /8tyM 3/2 , and v e is another con- 
stant called the velocity scale height, v e = (E/M) 1 ^ 2 ; both 
constants are determined by the total mass M and kinetic 
energy E of the ejecta, derived by integrating the mass den- 
sity and the kinetic energy density over space. A comparison 
with the power-law model yields the approximate power in- 
dex of the exponential density, n = —dlnp/dlnr = v/v e = 
r/v e t. Therefore, the exponential models with higher veloc- 
ity scale (increasing with E/M ratio) bear a flatter density 
distribution at a given flow velocity, and the density gradient 
increases with radius while it decreases over time. Such an 
exponential profile closely approximates the outer density 
distribution of the free ejecta that is obtained from many 
successful explosion models (Dwarkadas & Chevalier 1998), 
and the properties of the inferred ejecta clump in the most 
favorable explosion scenario, including the W7 deflagration 
model and the delayed detonation models, which depends 
on a near-Chandrasekhar-mass C-0 white dwarf explosion, 
are compatible with those needed by Tycho's knots (Wang 
2008). 

This study used an explosion mass of Mi. 4=1 and an 
explosion energy of Ebi=1, where Mi. 4 represents the mass 
in terms of the Chandrasekhar mass, 1.4Mq, and E51 de- 
notes the energy in units of 10 ergs, along with a constant 
ambient density of p am =2.34 x 10~ 24 gm cm -3 (correspond- 
ing to a hydrogen number density no = 1 cm -3 with a 
H/He ratio of 10/1 by number). Since the deceleration of 
an SNR is caused by the ambient medium, the interaction 
of the SN ejecta with a denser surrounding medium prompts 
a later evolutionary phase for the remnant. Nonetheless, so- 
lutions based on this set of parameters can be scaled into 
non-dimensional solutions, which can be re-scaled to ob- 
tain other dimensional solutions corresponding to various 
sets of M, E and p am - One evolutionary sequence in the 
non-dimensional variables thus represents virtually all pos- 
sible dimensional solutions. The following uses dimension- 
less variables, r' = r/R, v' = v/V and t' = t/T, where 
R' = (3A//4ti7w)5, V = (2£/M)5, and T" = R'/V are 
the scaling parameters, as described in WC01 to express 
the solution. For M1.4 = 1 and E51 = 1, R = 2.19 pc, 
V = 8.45 x 10 3 kms" 1 , and T = 248 yr. 
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Figure 1. Radial profile of gas density in one-dimensional expo- 
nential model at dimensionless age t' = 1.75 for 7=5/3, 4/3, 1.2, 
and 1.1, corresponding to Tycho's present state. Ratios of forward 
shock radius/reverse shock radius to reverse shock radius/forward 
shock radius are 1.51/0.66, 1.33/0.71, 1.22/0.81, and 1.15/0.87, 
respectively. 



4 EVOLUTION OF INSTABILITIES 

4.1 One-dimensional Simulations 

One-dimensional hydrodynamic simulations of the cases 
7=5/3, 4/3, 1.2, and 1.1 were performed first. Figure [l] 
presents the one-dimensional density profile in an age t' = 
1.75, which is approximately the normalized age of Tycho's 
remnant, to = 1-77 (438 years), if standard explosion param- 
eters are adopted. If a W7-like explosion condition, Mi. 4 = 1 
and E 5 i = 1.3, is assumed, then Tycho's true age corre- 
sponds to a dimensionless age of to = 2.02. Figures 2, 3, and 
4 illustrate the density distributions in the case of 7 = 4/3, 
1.2, and 1.1 at four various stages, overplotted with the 
angle-averaged radial profile. For 7 = 5/3, the reverse shock 
begins to turn inward in the stellar frame at t' = 2.5, and 
reaches the stellar center at t' = 8. In the non-adiabatic 
cases, because the sound speeds of the hot, shocked ejecta 
and hot, shocked ISM are still high and permit no pres- 
sure difference at their interface, the pressure-equilibrated 
contact discontinuity remains fixed in space while the inter- 
shock layer is further compressed toward the contact discon- 
tinuity. Accordingly, the turnover of the reverse shock and 
its subsequent impact with the explosion center are signifi- 
cantly delayed. X-ray observations have suggested that the 
reverse and forward shocks in Tycho's remnant are located 
at radii of 2.08 pc (r' = 0.98 in dimensionless units) and 3.09 
pc (r' = 1.45), and decelerate with deceleration parameters 
of 0.15 and 0.47, respectively. (See references in Dwarkadas 

6 Chevalier 1998 and Warren et al. 2005). Notably, these 
observed position and deceleration of the shocks are consis- 
tent with the result obtained using the exponential model 
with 7 = 5/3, but differ greatly from those obtained with 

7 = 1.1; at 7 = 5/3, the ratio of reverse to forward shock 
radii is ~ 2/3, whereas 7 = 1.1 yields an extremely high 
value of 0.87. The deviation of the shock radius ratios from 
the observed values increases with time because the use of 
small 7 causes the SNR to suffer from energy losses during 
expansion. 
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Figure 2. Angle-averaged two-dimensional density distribution 
plotted with the one-dimensional unperturbed solution at four 
stages of the evolution with 7 = 4/3. The time and radius use 
the normalized units given in the text. Two-dimensional solutions 
apply a 20% perturbation at the contact discontinuity. 
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Figure 3. Angle-averaged two-dimensional density distribution 
plotted with one-dimensional unperturbed solution at four stages 
of evolution with 7 = 1.2. Two-dimensional solutions apply a 20% 
perturbation below the reverse shock. 



4.2 Two-dimensional Simulations 

In most two-dimensional simulations, a density perturba- 
tion with an amplitude of 20% across contact discontinuity 
is applied to excite the instability as efficiently as possi- 
ble. The instability was seeded early enough for the flow 
to evolve into the nonlinear regime before Tycho's present 
epoch, when the growth of R-T fingers is matched by their 
destruction due to shearing and advection. Figures 5 and 6 
plot the density contours for 7 = 4/3 and 7 = 1.2, respec- 
tively. Following an initial growth phase, the initial pertur- 
bation grows to form mushroom-shaped caps. Shrinking of 
the intershock layer and increasing in the density contrast 
across the CD. generate narrower, denser and sharper R- 
T fingers. The flow then becomes independent of its initial 
conditions and enters a nonlinear phase, in which an increas- 
ingly larger instability mode is built up. The reverse shock 
front becomes corrugated and moves slightly inward, and 
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Figure 4. Angle-averaged two-dimensional density distribution 
plotted with the one-dimensional unperturbed solution at four 
stages of the evolution with 7 = 4/3. The time and radius use the 
normalized units given in the text. Two-dimensional solutions ap- 
ply a 20% perturbation at the contact discontinuity. The reverse 
and forward shocks are strongly perturbed and smear to an inner 
and outer radius respectively. 



the forward shock is not much perturbed. The evolution of 
the flow appears similar to that seen in WC01. In the case 
of 7 = 1.2, the R-T fingers can exert strong pressure and 
slightly disrupt the blast wave outline only in the initial 
stage of linear growth. However, the fingers at saturation 
do not come close to the blast wave, and so do not affect 
them; once the dense, shocked ejecta is mixed close to the 
forward shock, it quickly drops behind. The forward growth 
of the caps is impeded by the drag of the flow. The relative 
motion of flows between a finger and its surroundings bends 
the stems and disrupts the flow. The mushroom cap even- 
tually falls off to the side. The remaining filaments are then 
swept back to the CD. In the later stages, Kelvin-Helmholtz 
instability takes over, creating vortex rings in the less dense 
regions that are left by the original mushroom caps. The vor- 
tex rings gradually come to dominate the spikes, and only 
the stems remain recognizable. 

In the case of 7 = 1.1 (Fig. 7), the R-T fingers quickly 
extend to the forward shock in the linear stage of instability 
development, pushing small regions ahead of the average 
shock radius and distorting the outer blast wave. The R- 
T fingers exert strong pressure and disrupt the blast wave 
outline. The extent of the distortion seem to depend on the 
wavelength of the perturbation, and also the time since the 
initial perturbation is applied. However, as in the case of 7 = 
1.2, the vortex rings eventually come to dominate globally, 
and the bumpy structures on the remnant outline smooth 
down over time. 

The radius of the R-T fingers is compared to the forward 
shock radius. The extension of the R-T fingers is greatest in 
the stage of linear instability growth. At 7=1.1, the fingers 
pass the nominal shock radius. At 7 = 4/3 and 7 = 1.2 (Fig. 
[6] however, the R-T finger reaches only ~ 87% of the blast 
wave radius in Tycho's present epoch. The extent of mixing 
is virtually identical to that in the case of 7 = 5/3 considered 
by WC01. For 7 = 1.2 with efficient perturbations, the R-T 
fingers reach ~ 93% of the blast wave radius (Fig. [6| and 
slightly perturb the shock (Fig. [7f , but the mixture is much 
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less evident than the dense gas immediately beneath the 
forward shock front. 

The amplitude, wavelength and initial age of the per- 
turbation were varied. WC01 demonstrated that the insta- 
bility in the nonlinear regime is insensitive to the initial 
strength of the perturbation. Dwarkadas (2000) used the 
VH-1 code to indicate that the instability that arose from 
numerical noise remained in the linear regime, because vis- 
cosity and diffusion tend to dampen small-wavelength in- 
stability modes, which outpace all other modes; but with 
a 2% perturbation, his results were consistent with WC01. 
When the reverse shock accelerates back towards the cen- 
ter, distinct Rayleigh- Taylor fingers that protrude outwards 
from the reverse shock front are formed, when the intershock 
instability is not fully developed, whereas in the nonlinear 
case, the mode of such convection is not obvious (Fig. 8). 
During of the work on WC01, the flow could be induced into 
a nonlinear phase without an initial perturbation, as long as 
the simulation was begun soon enough with sufficiently high 
resolution, when the intershock structure had a high density 
contrast and large expansion rate. Given that Tycho's rem- 
nant shell has propagated for numerous inflight times since 
the explosion, the nonlinear instability should have become 
fully developed in the remnant. 

The use of 7 = 1.1 easily generates sufficiently strong 
excitations to influence the blast wave, particularly in the 
early stage of an SNR. The cases with 1% perturbation 
initiated at 10 4 s and with 20% perturbation initiated at 
10 6 s both yield strong turbulence, such that the bulges 
on the remnant outline change rapidly in space and time. 
The presence of the strong convection is not determined by 
the strength of the applied perturbation, but depends on 
whether there is a sufficiently large density contrast across 
the CD. Notably, because the fingers reach the forward 
shock before saturation, the strong shearing flow that is gen- 
erated while the blast wave is distorted blocks the growth 
of R-T fingers, which therefore cannot be pushed out be- 
yond the shock, and their growth stalls. As shown in BE01, 
if a fragment (R-T finger) of dense ejecta has sufficient ra- 
dial momentum to push out the forward shock in a local 
protrusion, then the deformed shock generates a substantial 
tangential post-shock shear flow around the ejecta finger, 
rapidly disrupting the finger or ejecta fragment through the 
Kelvin-Helmholtz instability, before it is advected back into 
the interaction region. Therefore, a local protrusion of the 
shock front of more than a few percent of the blast wave ra- 
dius cannot be formed when the unstable flow has saturated 
into the nonlinear regime. 



5 COMPARISON BETWEEN THE 
EXPONENTIAL MODEL AND THE 
POWER-LAW MODEL 

The exponential model can be understood as a power-law 
model with a continuously decreasing power index n. In 
BEOl's power-law model, the mixture with the highest com- 
pression ratio (<r = 21 and 7=1.1) quickly reaches and sur- 
passes the average shock radius; the remnant outline are 
perturbed with low-amplitude, small-wavelength bumps, al- 
though the average shock front remains almost spherical. 
Higher compression ratios is expected to cause significant de- 




Figure 6. Density images of the 7 = 1.2 model. Left panel - 
Simulation was initiated at t' = 0.0128 (1 X 10 s s). Initial density 
perturbation was imposed at 2.5 X 10 s s to ejecta below reverse 
shock, with ~ 20% amplitude and an I = 100 mode of the spher- 
ical harmonic function. Right panel - Simulation was initiated at 
t' = 0.00025 (1 X 10 6 s). Initial density perturbation was imposed 
at 4 X 10 6 s across contact discontinuity, with ~ 20% amplitude 
and an I = 100 mode of the spherical harmonic function. The grid 
has 600 uniform radial zones by 300 angular zones in one half of 
a quadrant. Contours correspond to logarithmically scaled values 
between the lowest and the highest values sampled in the inter- 
shock region. Color bars represent base-10 logarithmic values of 
density. 

viation from spherical symmetry. In the exponential model, 
a preferred mode in increasing wavelength that is indepen- 
dent of the initial conditions develops toward the end of the 
evolution. The perturbation to the remnant outline is even- 
tually dominated by the standard evolutionary paradigm 
that was illustrated in WC01 and therefore depends on the 
age of the remnant. A strong perturbation that is seeded 
near Tycho's age may protrude out of the blast wave, but 
this phenomenon is caused by the fact that the instabil- 
ity has not grown into the nonlinear regime. Given that 
a fully-grown perturbation causes less distortion than the 
perturbation in the linear stage, if the instability develops 
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Figure 5. Series of density contours that depict time evolution of dynamical instability using 7 = 3/4. Simulation was initiated at 
t' = 0.00025 (1 X 10 6 sec). Initial density perturbation was imposed at t' = 0.0006 (4 X 10 6 sec) to ejecta across the contact discontinuity, 
with 20% amplitude in the I = 100 mode of the spherical harmonic function. The grid has 600 nonuniform radial zones by 300 angular 
zones in one half of a quadrant, which resolved the intershock region into ~ 400 radial zones. Contours correspond to logarithmically 
scaled values between the lowest and the highest values sampled in the intershock region. Color bars represent base-10 logarithmic values 
of density. 







from a small perturbation soon after the SN explodes, then 
the remnant spherical outline should not be substantially 
altered. The result is consistent with that in the dramatic 
case of n — 12, s — (for a constant ambient medium, 
pam oc r~ a ), and 7 = 1.1, considered by BE01, in which 
the deviation from a spherical outer shock is small. Notably, 
the power index n of the free ejecta underneath the reverse 
shock drops below 7 at t' ^ 0.1, and the present density 
power index at Tycho's reverse shock is only n ~ 3 (see Fig. 
5 of WC01). In the power-law model self-similarity does not 
exist for n 3. 

BE01 showed that, despite the large changes in the ra- 
dial profiles, the width of the mixing region in terms of the 
ejecta mass fraction in the power-law model did not notice- 
ably vary with 7. The extent of mixing or the saturated 
instability growth rate, measured by summing the angular 



components of the turbulent energy density in the intershock 
region and then normalizing the sum to the bulk kinetic en- 
ergy density that is associated with the forward shock front 
(to remove the expansion and compression effects in the ra- 
dial direction), is also relatively unaffected by changes in 7. 
Hence, the dominant wavelength of the instability appears 
qualitatively similar for 7 = 4/3 and 7 = 1.2 (see Fig. 4 in 
BE01, for n = 7 and a circumstellar medium s — 2). In the 
exponential model, the maximum radius of the R-T fingers 
(in terms of the forward shock radius) does not vary sig- 
nificantly with 7. Furthermore, both 7 = 4/3 and 7 = 1.2 
are associated with a similar mode (wavenumber of ~ 4) to 
that in the adiabatic case of WC01 in Tycho's present epoch. 
The wavenumber is less than that in the adiabatic, power- 
law model with n — 7 and s = 0, as expected, where a half 
quadrant has six waves (Chevalier, Blondin, & Emmering, 
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Figure 8. Series of density contours that depict time evolution of dynamical instability using 7 = 1.1. Simulation was initiated at 
t' = 0.00025 (1 X 10 6 s). Initial density perturbation was imposed at t' = 0.0006 (4 X 10 6 s) to ejecta below reverse shock, with 20% 
amplitude in the / = 100 mode of the spherical harmonic function. The grid has 600 nonuniform radial zones by 300 angular zones in one 
half of a quadrant, which resolved the intershock region into ~ 400 radial zones. Contours correspond to logarithmically scaled values 
between the lowest and the highest values sampled in the intershock region. Color bars represent base-10 logarithmic values of density. 





1992). The instability has a delayed response to the change 
in the density gradient. The similarity among the dominant 
instability modes at various 7 above 1.2 follows from the fact 
that the energy losses in a SNR cannot change the radius 
of the contact discontinuity and so considerably reduce the 
decelerations of the shocks in Tycho's present epoch. This 
result is consistent with that of BE01. Notably, the case of 
7 = 1.2 yield a 12% lower speed of the blast wave than does 
7 = 5/3, whereas the case of 7 = 1.1 gives a 28% lower shock 
speed, amounting to unrealistic losses in kinetic energy. 

A major difference between the two models is that, in 
the power-law model, the instability is insensitive to changes 
in the compression ratio when the growth rate of R-T fin- 
gers (curves in Fig. 5 of BE01) become approximately flat, as 
when the system is in the nonlinear regime, while in the ex- 
ponential model, the instability decays continuously. Figure 



[To] plots the evolution of the turbulence, measured as inte- 
grating the angular component of the kinetic energy density 
over the volume of shocked gas; 



X = 



J pugdr 

JpdT 



(3) 



After a strong linear growth, the parameters converge to 
similarly small values of ~ 0.2 at t' = 2 and ~ 0.02 at 
t' — 8. Since the ratio of the forward shock radius to the 
contact discontinuity radius increases with time, the R-T 
fingers eventually drop behind the forward shock, and insta- 
bilities fade. If the power-law model is extended by including 
a flat component to the inner ejecta and if the system evolves 
long enough such that the reverse shock begins to turn over 
to the stellar center (Fig. 1 of WC02), then the ratio of 
the forward shock radius to the contact discontinuity radius 
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Figure 7. Density images of the 7 = 1.2 model at t' = 2.07, 
5.02, 9.12, and 11.69. Initial density perturbation was imposed at 
4 X 10 e s across contact discontinuity, with 20% amplitude and 
an I = 100 mode of the spherical harmonic function. The grid has 
600 uniform radial zones by 300 angular zones in one quadrant. 




Figure 9. Density images of the 7 = 1.1 models that depict 
different results at t' = 7.83, when the reverse shock has turned 
over. Contours correspond to logarithmically scaled values be- 
tween the lowest and the highest values sampled in the intershock 
region. Left - Perturbation was imposed with 10% amplitude and 
an I = 100 mode of the spherical harmonic function. The grid 
has 300 radial by 600 angular zones. Right - No perturbation was 
imposed. The grid has 600 radial by 100 angular zones. The in- 
stability at the reverse shock displays a dominant mode because 
the instability in the intershock region is not fully grown. Such a 
result is spurious. 



(which value remains fixed in the self-similar phase) starts 
to increase, and so the R-T convection declines as well. 



6 DISCUSSION AND CONCLUSIONS 

The acceleration by young SNR shocks of cosmic rays with 
high efficiency is expected to produce compression ratios of 
shocks that is considerably more than 4. To examine insta- 
bility in a real situation in which the shocked gas is cooled by 
the escape of superthermal cosmic-ray particles from the in- 
tershock region, low adiabatic indices are input to the expo- 
nential density model, which is applicable to Type la SNRs. 
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7 = 1.2 
7 = 4/3 
7 = 5/3 



1 T 



0.1 



0.01 



t/T 

Figure 10. Evolution of angular turbulent energy density nor- 
malized to bulk mass in intershock region for four values of 7. 
Same angular resolution of 300 zones in one half of a quadrant 
is used in all cases. Turbulence is measured in cgs units. Curve 
cross occurs when R-T fingers create bulges on blast wave. 



Since the extent of R-T mixing is determined by the de- 
celeration and the density profile of the intershock region, 
such a simple model should be able to represent the overall 
dynamics of the R-T instability that develops in a typical 
Type la SNR. 

When the adiabatic index is reduced, the SNR inter- 
shock gas is compressed toward the contact discontinuity, 
causing SN ejecta material to move closer to, or even slightly 
ahead of, the nominal shock radius. However, as the ratio of 
the forward shock radius to the contact discontinuity radius 
increases with time, the intershock structure becomes less 
compressed, and the R-T instabilities gradually decline. Pro- 
trusion from the shocks (which can be caused by a limited 
grid resolution) is expected only in the linear growth stage of 
the R-T instability, probably soon after the SN explodes. In 
later dynamic phases of a remnant, Kelvin-Helmholtz insta- 
bility overpowers Rayleigh- Taylor instability, and turbulence 
appearing in vortices dominates throughout the rest of the 
evolution, leaving the intershock flow rather featureless. 

When the shock compression reaches ~ 21 (7=1.1), the 
intershock region contracts to the extent that the reverse 
shock fronts resides at ^ 87% of the forward shock radius in 
the one-dimensional unperturbed case, or ~ 77% in the two- 
dimensional angle-averaged profile. However, a much smaller 
radius fraction of ~ 2/3 is observed for Tycho's SNR, which 
is consistent with adiabatic expansion. Contrary to previ- 
ous thought, the compression does not alter the degree of 
convective mixing; strong perturbations rapidly stall in the 
linear instability growth stage. For a remnant near Tycho's 
present epoch, the mixing takes place within a region below 
~ 87% of the blast wave radius, independently of the value 
of 7. This result substantiates BEOl's finding that the sat- 
urated instability growth rates in the power-law model do 
not detectably vary with 7. 

The changes in the dynamics of the convective flow that 
are caused by the energy losses for low values of 7 are small 
when 7 1.2. The extreme case of 7 = 1.1 involves sig- 
nificant, and probably not physically reasonable, losses of 
energy from the shock to cosmic rays. The escape of cosmic 
rays from shocks (by Alfven heating, among other processes) 
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is not well understood, and detailed models of cosmic-ray- 
modified shock remain theoretical. In Ferrand et al. (2010), 
the kinetic model of Blasi (2004), which solves the particle 
distribution and the fluid velocity as functions of the parti- 
cle momentum, is used to calculate the time-evolving shock 
properties and particle spectrum; a compact shock structure 
is revealed, but the R-T mixing remains not significantly af- 
fected. A magneto-hydrodynamic treatment of cosmic rays 
given in Skilling (1972 & 1975), which includes the parti- 
cles as a second fluid component in a plasma system, can 
be utilized for future modeling of cosmic-ray modified SNRs 
(for applying such governing equations to Parker instabili- 
ties please see Lo, Ko, & Wang 2011). 

This study supports our earlier determination that en- 
hanced R-T mixing is not responsible for the filaments and 
bumps seen on the outlines of young SNRs like Tycho's SNR. 
Notably, the above result also holds for core-collapse SNRs, 
as determined using the power-law model after t' 0.7, cor- 
responding to ~ 173 yr for Mi. 4 = 1 and E51 = 1, or ~ 890 
yr for a ten-fold explosion mass, when the reverse shock be- 
gins to enter the innermost plateau of the ejecta, terminat- 
ing the self-similar expansion phase (WC02). The conclusion 
is based on the dynamics of two-dimensional flows. Three- 
dimensional simulations in BE01 reveal more deformation of 
the forward shock, with an overall drop in the average com- 
pression ratio and an increase in the width of the intershock 
region, but the convection remains relatively unaffected by 
the dimensionality of the simulation. 

The contraction of the intershock structure depends on 
the pressure inside. External factors may adversely influence 
the flow. For example, if a strong magnetic field is present 
in the ISM, then the pressure from the shocked toroidal 
magnetic field will distend the forward-shocked layer, pre- 
venting the forward shocked layer from shrinking, and fur- 
ther impeding the R-T mixing. Radio polarization measure- 
ments of several supernova remnants however show the dom- 
inance of the radial magnetic field (Dickel, ven Breugel, & 
Milne 1991 for Tycho's SNR; Dickel, Strom, & Milne 2001 
for G315.4-23). These observations also reveal the presence 
of out-moving Mach cones that resembles the features in 
the Vela SNR (Strom et al. 1995; Aschenbach, Egger, & 
Truiimper 1995). On the other hand, if the forward shock 
encounters a clumpy interstellar medium, then the vortices 
that are generated behind the shock front can channel the R- 
T instability to enhance the mixing and create protrusions in 
the supernova remnant shell (Jun, Jones, & Norman 1996). 
Fujita et al. (2009) suggested that particles can be acceler- 
ated close to molecular clouds even after the SNR becomes 
radiative. However, the process for the Rayleigh- Taylor fin- 
gers reach the forward shock (depending on the properties 
and distribution of the clouds that are engulfed by the su- 
pernova shock) is transient, and whether the shock-cloud 
interaction can overcome the global decay of instability in 
a remnant like Tychos and in a typical, turbulent ISM is 
unclear. 

The possibility that the mixture of heavy-element ejecta 
near the blast wave is produced by R-T convection during 
the ejecta-dominated stage of an SNR is excluded, and we 
speculate that the injection of cosmic-ray particles at Ty- 
cho's forward shock is not sufficiently strong to alter the 
shock properties, including the compression ratio. The ex- 
treme case of 7 = 1.1 suggests that the outline is not to be 



substantially perturbed. Even if the injection and efficient 
acceleration of cosmic rays at the shock suffice, the extent of 
mixing is not expected to be significantly enhanced relative 
to the adiabatic case. 

While our result contradicts the interpretation of War- 
ren et al. (2005) based on X-ray data, the simulations simply 
reflect the basic dynamic properties of flow in multi dimen- 
sions that X-ray power spectra may not adequately reveal. 
The increased mixing is best explained by the interaction 
clumpy ejecta that is produced by the nickel bubble effect 
with the intershock region of the remnant. In this scenario, 
the radioactive heating of 56 Ni inflates the central ejecta 
and drives a strong thin shock to compress the outer heavy 
elements into dense ejecta clumps. Since the heavy elements 
were shocked before they underwent reverse shock, this pro- 
cess is likely to facilitate the creation of more energetic par- 
ticles, as well as denser clumping, in a supernova remnant. 
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